print("R codes for replicating ``Rejecting Non-Paternalist Motivation: An Experimental Test'' by Xianwen Chen and Øivind Schøyen")
print("Journal of Experimental Political Science")
print("Author of the Codes: Xianwen Chen")
print("Email: xianwen.chen@gmail.com")
print("Last revised on: April 1, 2021")

print("######################################")
print("######################################")

print("Load libraries.")

if (! require(dplyr) ) {
    install.packages("dplyr", repos = "https://cloud.r-project.org")
    library(dplyr)
}

if (! require(ggplot2) ) {
    install.packages("ggplot2", repos = "https://cloud.r-project.org")
    library(ggplot2)
}

if (! require(ks) ) {
    install.packages("ks", repos = "https://cloud.r-project.org")
    library(ks)
}

print("######################################")
print("######################################")

Spectator_Data <- read.csv( "01-data-spectators.csv" )
print("Data is read as the Spectator_Data data frame.")
print("Each row represents data from a spectator participant.")
print("The If_Treated variable took the value of 1 if the participant received treatment, and 0 otherwise.")
print("The If_Redist variable took the value of 1 if the participant preferred redistribution of income, and 0 otherwise.")
print("The Preva_Est variable recorded the participant's estimate of the prevalence of his or her redistribution preference.")
print("The If_WTP variable took the value of 1 if the participant was willing to implement his or her redistribution preference, and 0 otherwise.")
print("The Preva_True variable recorded the true prevalence of the participant's redistribution preference among the workers.")

print("######################################")
print("######################################")

print("Calculate the results presented in Table 1.")

print("Fairness views among workers were not calculated from the spectator data set, but from the worker data set. Thus the 33% and the 67% of fairness views are treated as given in this analysis.")

print("Total sample size:")
print(length(Spectator_Data[,1]))

print("Number of participants that did not receive information treatment and preferred redistribution:")
Spectator_Data_No_Treatment = Spectator_Data[Spectator_Data$If_Treated == 0, ]
Spectator_Data_No_Treatment_Yes_Distribution = Spectator_Data_No_Treatment[Spectator_Data_No_Treatment$If_Redist == 1, ]
print(length(Spectator_Data_No_Treatment_Yes_Distribution[,1]))

print("Among the participants that did not receive information treatment and preferred redistribution, how many were willing to pay to implement their redistribution preference:")
Spectator_Data_No_Treatment_Yes_Distribution %>%
    summarise_at(vars(If_WTP), funs(mean(., na.rm=FALSE))) %>%
    print()

print("Among the participants that did not receive information treatment and preferred redistribution, the average estimated prevalence of those who were not willing to pay (If_WTP == 0) and those who were willing to pay (If_WTP == 1):")
Spectator_Data_No_Treatment_Yes_Distribution %>%
    group_by(If_WTP) %>%
    summarise_at(vars(Preva_Est), funs(mean(., na.rm=FALSE))) %>%
    print()

print("Number of participants that did not receive information treatment and preferred no redistribution:")
Spectator_Data_No_Treatment_No_Distribution = Spectator_Data_No_Treatment[Spectator_Data_No_Treatment$If_Redist == 0, ]
print(length(Spectator_Data_No_Treatment_No_Distribution[,1]))

print("Among the participants that did not receive information treatment and preferred no redistribution, how many were willing to pay to implement their redistribution preference:")
Spectator_Data_No_Treatment_No_Distribution %>%
    summarise_at(vars(If_WTP), funs(mean(., na.rm=FALSE))) %>%
    print()

print("Among the participants that did not receive information treatment and preferred no redistribution, the average estimated prevalence of those who were not willing to pay (If_WTP == 0) and those who were willing to pay (If_WTP == 1):")
Spectator_Data_No_Treatment_No_Distribution %>%
    group_by(If_WTP) %>%
    summarise_at(vars(Preva_Est), funs(mean(., na.rm=FALSE))) %>%
    print()

print("Number of participants that received information treatment and preferred redistribution:")
Spectator_Data_Yes_Treatment = Spectator_Data[Spectator_Data$If_Treated == 1, ]
Spectator_Data_Yes_Treatment_Yes_Distribution = Spectator_Data_Yes_Treatment[Spectator_Data_Yes_Treatment$If_Redist == 1, ]
print(length(Spectator_Data_Yes_Treatment_Yes_Distribution[,1]))

print("Among the participants that received information treatment and preferred redistribution, how many were willing to pay to implement their redistribution preference:")
Spectator_Data_Yes_Treatment_Yes_Distribution %>%
    summarise_at(vars(If_WTP), funs(mean(., na.rm=FALSE))) %>%
    print()

print("Among the participants that received information treatment and preferred redistribution, the average estimated prevalence of those who were not willing to pay (If_WTP == 0) and those who were willing to pay (If_WTP == 1):")
Spectator_Data_Yes_Treatment_Yes_Distribution %>%
    group_by(If_WTP) %>%
    summarise_at(vars(Preva_Est), funs(mean(., na.rm=FALSE))) %>%
    print()

print("Number of participants that received information treatment and preferred no redistribution:")
Spectator_Data_Yes_Treatment_No_Distribution = Spectator_Data_Yes_Treatment[Spectator_Data_Yes_Treatment$If_Redist == 0, ]
print(length(Spectator_Data_Yes_Treatment_No_Distribution[,1]))

print("Among the participants that received information treatment and preferred no redistribution, how many were willing to pay to implement their redistribution preference:")
Spectator_Data_Yes_Treatment_No_Distribution %>%
    summarise_at(vars(If_WTP), funs(mean(., na.rm=FALSE))) %>%
    print()

print("Among the participants that received information treatment and preferred no redistribution, the average estimated prevalence of those who were not willing to pay (If_WTP == 0) and those who were willing to pay (If_WTP == 1):")
Spectator_Data_Yes_Treatment_No_Distribution %>%
    group_by(If_WTP) %>%
    summarise_at(vars(Preva_Est), funs(mean(., na.rm=FALSE))) %>%
    print()

print("######################################")
print("######################################")

print("Calculate the results presented in Table 2.")

print("Model (1)")
Spectator_Data$If_Over_Projection <- 0
Spectator_Data$If_Over_Projection[ Spectator_Data$Preva_Est > Spectator_Data$Preva_True] <- 1
Reg_WTP_Overly_Projected <- lm( formula = If_WTP ~ If_Over_Projection, data = Spectator_Data )
print(summary( Reg_WTP_Overly_Projected ))

print("Model (2)")
Reg_WTP_Treatment <- lm( formula = If_WTP ~ If_Treated, data = Spectator_Data )
print(summary( Reg_WTP_Treatment ))

print("Model (3)")
Spectator_Data$If_Over_Projection_and_If_Treated <- Spectator_Data$If_Over_Projection * Spectator_Data$If_Treated
Reg_WTP_Treatment_Interaction_Overly_Projected <- lm( formula = If_WTP ~ If_Over_Projection + If_Treated + If_Over_Projection_and_If_Treated, data = Spectator_Data )
print(summary( Reg_WTP_Treatment_Interaction_Overly_Projected ))

print("######################################")
print("######################################")

print("Calculate the results presented in Table 3 in Appendix")

print("Pooled sample")
Reg_WTP_Correlation_with_WTP <- lm( formula = If_WTP ~ Preva_Est, data = Spectator_Data )
print(summary( Reg_WTP_Correlation_with_WTP ))

print("Control group")
Reg_WTP_Correlation_with_WTP_Control_Group <- lm( formula = If_WTP ~ Preva_Est, data = Spectator_Data_No_Treatment )
print(summary( Reg_WTP_Correlation_with_WTP_Control_Group ))

print("Treatment group")
Reg_WTP_Correlation_with_WTP_Treatment_Group <- lm( formula = If_WTP ~ Preva_Est, data = Spectator_Data_Yes_Treatment )
print(summary( Reg_WTP_Correlation_with_WTP_Treatment_Group ))

print("######################################")
print("######################################")

print("Calculate the results presented in Table 4 in Appendix")

print("All spectators")
Reg_WTP_Treatment_Effect <- lm( formula = If_WTP ~ If_Treated, data = Spectator_Data )
print(summary( Reg_WTP_Treatment_Effect ))

print("Those preferred redistribution")
Spectator_Data_Yes_Redistribution = Spectator_Data[Spectator_Data$If_Redist == 1, ]
Reg_WTP_Treatment_Effect_Yes_Redistribution <- lm( formula = If_WTP ~ If_Treated, data = Spectator_Data_Yes_Redistribution )
print(summary( Reg_WTP_Treatment_Effect_Yes_Redistribution ))

print("Those preferred no redistribution")
Spectator_Data_No_Redistribution = Spectator_Data[Spectator_Data$If_Redist == 0, ]
Reg_WTP_Treatment_Effect_No_Redistribution <- lm( formula = If_WTP ~ If_Treated, data = Spectator_Data_No_Redistribution )
print(summary( Reg_WTP_Treatment_Effect_No_Redistribution ))

print("######################################")
print("######################################")

print("Calculate the results presented in Table 5 in Appendix")

print("Projection bias > 20%")
Spectator_Data$Preva_Diff <- Spectator_Data$Preva_Est - Spectator_Data$Preva_True
Spectator_Data_Large_Prev_Diff_20 <- Spectator_Data[ Spectator_Data$Preva_Diff > 20, ]
Reg_WTP_Correlation_with_WTP_Large_Prev_Diff_20 <- lm( formula = If_WTP ~ Preva_Est, data = Spectator_Data_Large_Prev_Diff_20 )
print(summary( Reg_WTP_Correlation_with_WTP_Large_Prev_Diff_20 ))

print("Projection bias > 30%")
Spectator_Data$Preva_Diff <- Spectator_Data$Preva_Est - Spectator_Data$Preva_True
Spectator_Data_Large_Prev_Diff_30 <- Spectator_Data[ Spectator_Data$Preva_Diff > 30, ]
Reg_WTP_Correlation_with_WTP_Large_Prev_Diff_30 <- lm( formula = If_WTP ~ Preva_Est, data = Spectator_Data_Large_Prev_Diff_30 )
print(summary( Reg_WTP_Correlation_with_WTP_Large_Prev_Diff_30 ))

print("Projection bias > 40%")
Spectator_Data$Preva_Diff <- Spectator_Data$Preva_Est - Spectator_Data$Preva_True
Spectator_Data_Large_Prev_Diff_40 <- Spectator_Data[ Spectator_Data$Preva_Diff > 40, ]
Reg_WTP_Correlation_with_WTP_Large_Prev_Diff_40 <- lm( formula = If_WTP ~ Preva_Est, data = Spectator_Data_Large_Prev_Diff_40 )
print(summary( Reg_WTP_Correlation_with_WTP_Large_Prev_Diff_40 ))

print("######################################")
print("######################################")

print("Calculate the results presented in Table 6 in Appendix")

print("Test (1)")
print(str(ks.test( x = Spectator_Data_No_Treatment$Preva_Est, y = Spectator_Data_Yes_Treatment$Preva_Est, alternative = "two.sided" )))
print("Sample sizes:")
print(length(Spectator_Data_No_Treatment$Preva_Est))
print(length(Spectator_Data_Yes_Treatment$Preva_Est))
print("")

print("Test (2)")
Spectator_Data_No_Treatment_Yes_WTP = Spectator_Data_No_Treatment[Spectator_Data_No_Treatment$If_WTP == 1, ]
Spectator_Data_Yes_Treatment_Yes_WTP = Spectator_Data_Yes_Treatment[Spectator_Data_Yes_Treatment$If_WTP == 1, ]
print(str(ks.test( x = Spectator_Data_No_Treatment_Yes_WTP$Preva_Est, y = Spectator_Data_Yes_Treatment_Yes_WTP$Preva_Est, alternative = "two.sided" )))
print("Sample sizes:")
print(length(Spectator_Data_No_Treatment_Yes_WTP$Preva_Est))
print(length(Spectator_Data_Yes_Treatment_Yes_WTP$Preva_Est))
print("")

print("Test (3)")
Spectator_Data_No_Treatment_Yes_WTP_Yes_Redistribution = Spectator_Data_No_Treatment_Yes_WTP[Spectator_Data_No_Treatment_Yes_WTP$If_Redist == 1, ]
Spectator_Data_Yes_Treatment_Yes_WTP_Yes_Redistribution = Spectator_Data_Yes_Treatment_Yes_WTP[Spectator_Data_Yes_Treatment_Yes_WTP$If_Redist == 1, ]
print(str(ks.test( x = Spectator_Data_No_Treatment_Yes_WTP_Yes_Redistribution$Preva_Est, y = Spectator_Data_Yes_Treatment_Yes_WTP_Yes_Redistribution$Preva_Est, alternative = "two.sided" )))
print("Sample sizes:")
print(length(Spectator_Data_No_Treatment_Yes_WTP_Yes_Redistribution$Preva_Est))
print(length(Spectator_Data_Yes_Treatment_Yes_WTP_Yes_Redistribution$Preva_Est))
print("")

print("Test (4)")
Spectator_Data_No_Treatment_Yes_WTP_No_Redistribution = Spectator_Data_No_Treatment_Yes_WTP[Spectator_Data_No_Treatment_Yes_WTP$If_Redist == 0, ]
Spectator_Data_Yes_Treatment_Yes_WTP_No_Redistribution = Spectator_Data_Yes_Treatment_Yes_WTP[Spectator_Data_Yes_Treatment_Yes_WTP$If_Redist == 0, ]
print(str(ks.test( x = Spectator_Data_No_Treatment_Yes_WTP_No_Redistribution$Preva_Est, y = Spectator_Data_Yes_Treatment_Yes_WTP_No_Redistribution$Preva_Est, alternative = "two.sided" )))
print("Sample sizes:")
print(length(Spectator_Data_No_Treatment_Yes_WTP_No_Redistribution$Preva_Est))
print(length(Spectator_Data_Yes_Treatment_Yes_WTP_No_Redistribution$Preva_Est))
print("")

print("Test (5)")
Spectator_Data_Yes_Treatment_No_WTP = Spectator_Data_Yes_Treatment[Spectator_Data_Yes_Treatment$If_WTP == 0, ]
print(str(ks.test( x = Spectator_Data_Yes_Treatment_Yes_WTP$Preva_Est, y = Spectator_Data_Yes_Treatment_No_WTP$Preva_Est, alternative = "two.sided" )))
print("Sample sizes:")
print(length(Spectator_Data_Yes_Treatment_Yes_WTP$Preva_Est))
print(length(Spectator_Data_Yes_Treatment_No_WTP$Preva_Est))
print("")

print("Test (6)")
Spectator_Data_Yes_Treatment_No_WTP_Yes_Redistribution = Spectator_Data_Yes_Treatment_No_WTP[Spectator_Data_Yes_Treatment_No_WTP$If_Redist == 1, ]
print(str(ks.test( x = Spectator_Data_Yes_Treatment_Yes_WTP_Yes_Redistribution$Preva_Est, y = Spectator_Data_Yes_Treatment_No_WTP_Yes_Redistribution$Preva_Est, alternative = "two.sided" )))
print("Sample sizes:")
print(length(Spectator_Data_Yes_Treatment_Yes_WTP_Yes_Redistribution$Preva_Est))
print(length(Spectator_Data_Yes_Treatment_No_WTP_Yes_Redistribution$Preva_Est))
print("")

print("Test (7)")
Spectator_Data_Yes_Treatment_No_WTP_No_Redistribution = Spectator_Data_Yes_Treatment_No_WTP[Spectator_Data_Yes_Treatment_No_WTP$If_Redist == 0, ]
print(str(ks.test( x = Spectator_Data_Yes_Treatment_Yes_WTP_No_Redistribution$Preva_Est, y = Spectator_Data_Yes_Treatment_No_WTP_No_Redistribution$Preva_Est, alternative = "two.sided" )))
print("Sample sizes:")
print(length(Spectator_Data_Yes_Treatment_Yes_WTP_No_Redistribution$Preva_Est))
print(length(Spectator_Data_Yes_Treatment_No_WTP_No_Redistribution$Preva_Est))

print("######################################")
print("######################################")

print("Create the four figures is Figure 1")

print("Calculate the histogram bins for those in the control group that preferred no-redistribution.")

Projected_Prevalence_Bin_No_Redistribution_Control_Group = data.frame( Count = c( 1:10 ), Bins =  c( "0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81-90", "91-100" ), Count_WTP = c( 1:10 ) )
Projected_Prevalence_Bin_No_Redistribution_Control_Group
Spectator_Data_No_Redistribution = Spectator_Data[ Spectator_Data$If_Redist == 0, ]
Spectator_Data_No_Redistribution_Control_Group = Spectator_Data_No_Redistribution[ Spectator_Data_No_Redistribution$If_Treated == 0, ]

i = 1
while ( i <= 10 ) 
{
	if ( i == 1 ) {
		Spectator_Data_Temp = Spectator_Data_No_Redistribution_Control_Group[ as.numeric( Spectator_Data_No_Redistribution_Control_Group$Preva_Est ) <= 10, ]
	} else {
		Spectator_Data_Temp = Spectator_Data_No_Redistribution_Control_Group[ as.numeric( Spectator_Data_No_Redistribution_Control_Group$Preva_Est ) >  ( ( i - 1)  * 10 ), ]
		Spectator_Data_Temp = Spectator_Data_Temp[ as.numeric( Spectator_Data_Temp$Preva_Est ) <= ( i * 10 ), ]
	}
	Projected_Prevalence_Bin_No_Redistribution_Control_Group$Count[ i ] = length(  Spectator_Data_Temp[ , 1] )
	Projected_Prevalence_Bin_No_Redistribution_Control_Group$Count_WTP[ i] = sum( Spectator_Data_Temp$If_WTP )
	i = i + 1
}

Projected_Prevalence_Bin_No_Redistribution_Control_Group
sum( Projected_Prevalence_Bin_No_Redistribution_Control_Group[ , 1 ] )
Projected_Prevalence_Bin_No_Redistribution_Control_Group$Percentage = Projected_Prevalence_Bin_No_Redistribution_Control_Group$Count / sum( Projected_Prevalence_Bin_No_Redistribution_Control_Group[ , 1 ] ) * 100
Projected_Prevalence_Bin_No_Redistribution_Control_Group

print("Create the histogram for those in the control group that preferred no-redistribution.")

Projected_Prevalence_Histograph_Manual_No_Redistribution_Control_Group = ggplot( Projected_Prevalence_Bin_No_Redistribution_Control_Group, aes( Bins, Count)) +  
	geom_col( fill = rgb( red = 255, green = 0, blue = 255, max = 255 ) ) 	+ 
	ylim(0, 200) + xlab("Estimated Prevalence of Fairness View (%)") + ylab("Number of Participants") + theme(axis.title.y = element_text(size = rel(1.8), angle = 90)) + theme(axis.title.x = element_text(size = rel(1.8), angle = 00)) + labs(title = "Find No-Redistribution Fair (Control Group)") + geom_col( aes( Bins, Count_WTP), fill = "blue" ) + theme( plot.title = element_text( size = 20 ) ) + theme(axis.title.y = element_text(size = rel(1.8), angle = 90)) + theme(axis.title.x = element_text(size = rel(1.8), angle = 00)) +
	geom_vline( xintercept = 7.2, linetype="dotted" ) +
	theme( axis.text=element_text(size=13) )
Projected_Prevalence_Histograph_Manual_No_Redistribution_Control_Group

print("The figure is saved as 03-output_from_R_codes-Histogram_Control_Group_No_Redistribution.pdf.")

ggsave("03-output_from_R_codes-Histogram_Control_Group_No_Redistribution.pdf")

print("Calculate the histogram bins for those in the control group that preferred redistribution.")

Projected_Prevalence_Bin_Yes_Redistribution_Control_Group = data.frame( Count = c( 1:10 ), Bins =  c( "0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81-90", "91-100" ), Count_WTP = c( 1:10 ) )
Projected_Prevalence_Bin_Yes_Redistribution_Control_Group
Spectator_Data_Yes_Redistribution = Spectator_Data[ Spectator_Data$If_Redist == 1, ]
Spectator_Data_Yes_Redistribution_Control_Group = Spectator_Data_Yes_Redistribution[ Spectator_Data_Yes_Redistribution$If_Treated == 0, ]

i = 1
while ( i <= 10 ) 
{
	if ( i == 1 ) {
		Spectator_Data_Temp = Spectator_Data_Yes_Redistribution_Control_Group[ as.numeric( Spectator_Data_Yes_Redistribution_Control_Group$Preva_Est ) <= 10, ]
	} else {
		Spectator_Data_Temp = Spectator_Data_Yes_Redistribution_Control_Group[ as.numeric( Spectator_Data_Yes_Redistribution_Control_Group$Preva_Est ) >  ( ( i - 1)  * 10 ), ]
		Spectator_Data_Temp = Spectator_Data_Temp[ as.numeric( Spectator_Data_Temp$Preva_Est ) <= ( i * 10 ), ]
	}
	Projected_Prevalence_Bin_Yes_Redistribution_Control_Group$Count[ i ] = length(  Spectator_Data_Temp[ , 1] )
	Projected_Prevalence_Bin_Yes_Redistribution_Control_Group$Count_WTP[ i] = sum( Spectator_Data_Temp$If_WTP )
	i = i + 1
}

Projected_Prevalence_Bin_Yes_Redistribution_Control_Group
sum( Projected_Prevalence_Bin_Yes_Redistribution_Control_Group[ , 1 ] )
Projected_Prevalence_Bin_Yes_Redistribution_Control_Group$Percentage = Projected_Prevalence_Bin_Yes_Redistribution_Control_Group$Count / sum( Projected_Prevalence_Bin_Yes_Redistribution_Control_Group[ , 1 ] ) * 100
Projected_Prevalence_Bin_Yes_Redistribution_Control_Group

print("Create the histogram for those in the control group that preferred redistribution.")

Projected_Prevalence_Histograph_Manual_Yes_Redistribution_Control_Group = ggplot( Projected_Prevalence_Bin_Yes_Redistribution_Control_Group, aes( Bins, Count)) +  
	geom_col( fill = rgb( red = 255, green = 0, blue = 255, max = 255 ) ) 	+ 
	ylim(0, 200) + xlab("Estimated Prevalence of Fairness View (%)") + ylab("Number of Participants") + theme(axis.title.y = element_text(size = rel(1.8), angle = 90)) + theme(axis.title.x = element_text(size = rel(1.8), angle = 00)) + labs(title = "Find Redistribution Fair (Control Group)") + geom_col( aes( Bins, Count_WTP), fill = "blue" ) + theme( plot.title = element_text( size = 20 ) ) +
	geom_vline( xintercept = 3.8, linetype="dotted" ) +
	theme( axis.text=element_text(size=13) )
Projected_Prevalence_Histograph_Manual_Yes_Redistribution_Control_Group

print("The figure is saved as 03-output_from_R_codes-Histogram_Control_Group_Yes_Redistribution.pdf.")

ggsave("03-output_from_R_codes-Histogram_Control_Group_Yes_Redistribution.pdf")

print("Calculate the histogram bins for those in the treatment group that preferred no-redistribution.")

Projected_Prevalence_Bin_No_Redistribution_Treatment_Group = data.frame( Count = c( 1:10 ), Bins =  c( "0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81-90", "91-100" ), Count_WTP = c( 1:10 ) )
Projected_Prevalence_Bin_No_Redistribution_Treatment_Group
Spectator_Data_No_Redistribution = Spectator_Data[ Spectator_Data$If_Redist == 0, ]
Spectator_Data_No_Redistribution_Treatment_Group = Spectator_Data_No_Redistribution[ Spectator_Data_No_Redistribution$If_Treated == 1, ]

i = 1
while ( i <= 10 ) 
{
	if ( i == 1 ) {
		Spectator_Data_Temp = Spectator_Data_No_Redistribution_Treatment_Group[ as.numeric( Spectator_Data_No_Redistribution_Treatment_Group$Preva_Est ) <= 10, ]
	} else {
		Spectator_Data_Temp = Spectator_Data_No_Redistribution_Treatment_Group[ as.numeric( Spectator_Data_No_Redistribution_Treatment_Group$Preva_Est ) >  ( ( i - 1)  * 10 ), ]
		Spectator_Data_Temp = Spectator_Data_Temp[ as.numeric( Spectator_Data_Temp$Preva_Est ) <= ( i * 10 ), ]
	}
	Projected_Prevalence_Bin_No_Redistribution_Treatment_Group$Count[ i ] = length(  Spectator_Data_Temp[ , 1] )
	Projected_Prevalence_Bin_No_Redistribution_Treatment_Group$Count_WTP[ i] = sum( Spectator_Data_Temp$If_WTP )
	i = i + 1
}

Projected_Prevalence_Bin_No_Redistribution_Treatment_Group
sum( Projected_Prevalence_Bin_No_Redistribution_Treatment_Group[ , 1 ] )
Projected_Prevalence_Bin_No_Redistribution_Treatment_Group$Percentage = Projected_Prevalence_Bin_No_Redistribution_Treatment_Group$Count / sum( Projected_Prevalence_Bin_No_Redistribution_Treatment_Group[ , 1 ] ) * 100
Projected_Prevalence_Bin_No_Redistribution_Treatment_Group

print("Create the histogram for those in the treatment group that preferred no-redistribution.")

Projected_Prevalence_Histograph_Manual_No_Redistribution_Treatment_Group = ggplot( Projected_Prevalence_Bin_No_Redistribution_Treatment_Group, aes( Bins, Count)) +  
	geom_col( fill = rgb( red = 255, green = 0, blue = 255, max = 255 ) ) 	+ 
	ylim(0, 200) + xlab("Estimated Prevalence of Fairness View (%)") + ylab("Number of Participants") + theme(axis.title.y = element_text(size = rel(1.8), angle = 90)) + theme(axis.title.x = element_text(size = rel(1.8), angle = 00)) + labs(title = "Find No-Redistribution Fair (Treatment Group)") + geom_col( aes( Bins, Count_WTP), fill = "blue" ) + theme( plot.title = element_text( size = 20 ) ) + theme(axis.title.y = element_text(size = rel(1.8), angle = 90)) + theme(axis.title.x = element_text(size = rel(1.8), angle = 00)) +
	geom_vline( xintercept = 7.2, linetype="dotted" ) +
	theme( axis.text=element_text(size=13) )
Projected_Prevalence_Histograph_Manual_No_Redistribution_Treatment_Group

print("The figure is saved as 03-output_from_R_codes-Histogram_Treatment_Group_No_Redistribution.pdf.")

ggsave("03-output_from_R_codes-Histogram_Treatment_Group_No_Redistribution.pdf")

print("Calculate the histogram bins for those in the treatment group that preferred redistribution.")

Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group = data.frame( Count = c( 1:10 ), Bins = c( "0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81-90", "91-100" ), Count_WTP = c( 1:10 ) )
Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group
Spectator_Data_Yes_Redistribution = Spectator_Data[ Spectator_Data$If_Redist == 1, ]
Spectator_Data_Yes_Redistribution_Treatment_Group = Spectator_Data_Yes_Redistribution[ Spectator_Data_Yes_Redistribution$If_Treated == 1, ]

i = 1
while ( i <= 10 ) 
{
	if ( i == 1 ) {
		Spectator_Data_Temp = Spectator_Data_Yes_Redistribution_Treatment_Group[ as.numeric( Spectator_Data_Yes_Redistribution_Treatment_Group$Preva_Est ) <= 10, ]
	} else {
		Spectator_Data_Temp = Spectator_Data_Yes_Redistribution_Treatment_Group[ as.numeric( Spectator_Data_Yes_Redistribution_Treatment_Group$Preva_Est ) >  ( ( i - 1)  * 10 ), ]
		Spectator_Data_Temp = Spectator_Data_Temp[ as.numeric( Spectator_Data_Temp$Preva_Est ) <= ( i * 10 ), ]
	}
	Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group$Count[ i ] = length(  Spectator_Data_Temp[ , 1] )
	Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group$Count_WTP[ i] = sum( Spectator_Data_Temp$If_WTP )
	i = i + 1
}

Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group
sum( Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group[ , 1 ] )
Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group$Percentage = Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group$Count / sum( Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group[ , 1 ] ) * 100
Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group

print("Create the histogram for those in the treatment group that preferred redistribution.")

Projected_Prevalence_Histograph_Manual_Yes_Redistribution_Treatment_Group = 
	ggplot( Projected_Prevalence_Bin_Yes_Redistribution_Treatment_Group, aes( Bins, Count)) +  
	geom_col( fill = rgb( red = 255, green = 0, blue = 255, max = 255 ) ) 	+ 
	ylim(0, 200) + 
	xlab("Estimated Prevalence of Fairness View (%)") + 
	ylab("Number of Participants") +
	theme(axis.title.y = element_text(size = rel(1.8), angle = 90)) + 
	theme(axis.title.x = element_text(size = rel(1.8), angle = 00)) + 
	labs(title = "Find Redistribution Fair (Treatment Group)") + 
	geom_col( aes( Bins, Count_WTP), fill = "blue" ) + 
	theme( plot.title = element_text( size = 20 ) ) +
	geom_vline( xintercept = 3.8, linetype="dotted" ) +
	theme( axis.text=element_text(size=13) ) 
Projected_Prevalence_Histograph_Manual_Yes_Redistribution_Treatment_Group

print("The figure is saved as 03-output_from_R_codes-Histogram_Treatment_Group_Yes_Redistribution.pdf.")

ggsave("03-output_from_R_codes-Histogram_Treatment_Group_Yes_Redistribution.pdf")


